#---- Load libraries and funtions -----

library(haven)
library(scales)
library(lme4)
library(lmerTest)
library(effects)
library(optimx)
library(emmeans)
library(multcomp)
library(multcompView)
library(ggeffects)
library(sjPlot)
library(sjstats)
library(sjmisc)
library(ggplot2)
library(gridExtra)
library(dplyr)
library(jtools)
library(car)
library(psych)
library(interactions)
library(lattice)
library(arm)
library(arsenal)
library(psych)

options(scipen=999)# non-scientific notation

# scale between 0-1#-----------------------------------

scale01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))} 

#reading data#--------------------------------------------

eu28 <- as.data.frame(read_sav('europe_covid_climate_2021.sav'))

#scaling and centering variables#--------------------------------------------
eu28$gender = scale01(eu28$gender)
eu28$gender_c = scale(eu28$gender,center=TRUE,scale=FALSE)
eu28$vacc_c = scale(eu28$vacc,center=TRUE,scale=FALSE)
eu28$age = scale01(eu28$age)
eu28$age_c = scale(eu28$age,center=TRUE,scale=FALSE)
eu28$isced = scale01(eu28$edu)
eu28$isced_c = scale(eu28$isced,center=TRUE,scale=FALSE)
describe(eu28$vacc_c)
eu28$leftright = scale01(eu28$leftright)
eu28$leftright_c = scale(eu28$leftright,center=TRUE,scale=FALSE)
eu28$fcv = scale01(eu28$fcv)
eu28$fcv_c = scale(eu28$fcv,center=TRUE,scale=FALSE)
eu28$settlement = scale01(eu28$settlement)
eu28$settlement_c = scale(eu28$settlement,center=TRUE,scale=FALSE)
eu28$hdi = scale01(eu28$hdi)
eu28$hdi_c = scale(eu28$hdi,center=TRUE,scale=FALSE)
eu28$cri = scale01(eu28$cri)
eu28$cri_c = scale(eu28$cri,center=TRUE,scale=FALSE)
eu28$totalcase = scale01(eu28$totalcase)
eu28$totalcase_c = scale(eu28$totalcase,center=TRUE,scale=FALSE)
eu28$totaldeath = scale01(eu28$totaldeath)
eu28$totaldeath_c = scale(eu28$totaldeath,center=TRUE,scale=FALSE)
eu28$vaccin = scale01(eu28$vaccin)
eu28$vaccine_c = scale(eu28$vaccin,center=TRUE,scale=FALSE)
eu28$str = scale01(eu28$stringency)
eu28$str_c = scale(eu28$str,center=TRUE,scale=FALSE)
eu28$cntry = as.factor(eu28$cntry)


#climate change concerns#--------------------------------------------

#null-model#--------------------------------------------

cwmodel0<-lmer(climateconcerns ~ 1+(1|cntry),data=eu28,REML=FALSE, weights = analysis_weight)

summ(cwmodel0)

#model1 - fixed-effects#-------------------------------------------


cwmodel1 = lmer(climateconcerns ~ 
                  fcv_c+vacc+leftright_c+gender + age_c + isced_c +
                  settlement_c+
                  totalcase_c+totaldeath_c+vaccine_c+str_c+
                  hdi_c+cri_c+
                  (1|cntry)
                , data = eu28, REML = FALSE,
                control=lmerControl(optimizer="optimx",
                                    optCtrl=list(method='nlminb')), 
                na.action = na.exclude,
                weights = analysis_weight)

summ(cwmodel1)


#model2 - interactions#-------------------------------------------

cwmodel2 = lmer(climateconcerns ~ 
                  fcv_c+vacc+leftright_c+gender + age_c + isced_c +
                  settlement_c+
                  totalcase_c+totaldeath_c+vaccine_c+str_c+
                  hdi_c+cri_c+fcv_c*leftright_c+fcv_c*age_c+
                  fcv_c*settlement_c+fcv_c*totalcase_c+fcv_c*str_c+
                  (1|cntry)
                , data = eu28, REML = FALSE,
                control=lmerControl(optimizer="optimx",
                                    optCtrl=list(method='nlminb')), 
                weights = analysis_weight)

summary(cwmodel2)

#model3 - random slope#-------------------------------------------

cwmodel3 = lmer(climateconcerns ~ 
                  fcv_c+vacc+leftright_c+gender + age_c + isced_c +
                  settlement_c+
                       totalcase_c+totaldeath_c+vaccine_c+str_c+
                       hdi_c+cri_c+fcv_c*leftright_c+fcv_c*age_c+
                  fcv_c*settlement_c+fcv_c*totalcase_c+fcv_c*str_c+
                       (1+fcv_c|cntry)
                     , data = eu28, REML = FALSE,
                     control=lmerControl(optimizer="optimx",
                                         optCtrl=list(method='nlminb')),
                na.action = na.exclude,
                     weights = analysis_weight)

summary(cwmodel3)

#computing likelihood ratio tests#--------------------------------------------

anova(cwmodel1, cwmodel2, refit = FALSE)
anova(cwmodel2, cwmodel3, refit = FALSE)


#plotting interactions#-------------------------------------------

describe(eu28$leftright_c)


plot1 <- interact_plot(model = cwmodel3, pred  = fcv_c,
                       modx  = leftright_c,
                       x.label = "FCV-19S",
                       modx.values = c(-0.4, 0.4),
                       modx.labels = c("Left","Right"),
                       y.label = "Climate change concerns",
                       legend.main = "Pol. orientation",
                       interval = TRUE)
plot11 <- (plot1 + theme(legend.position = "top",
                         panel.grid.major = element_blank(),
                         panel.grid.minor = element_blank(),
                         legend.title = element_text (size =9), 
                         legend.text = element_text (size=9), 
                         legend.key.width=unit(0.8,"cm")))

plot11

describe(eu28$age_c)

plot2 <- interact_plot(model = cwmodel3, pred  = fcv_c,
                       modx  = age_c,
                       x.label = "FCV-19S",
                       modx.values = c(-0.2, 0.2),
                       modx.labels = c("Younger","Older"),
                       y.label = "",
                       legend.main = "Age",
                       interval = TRUE)
plot22 <- (plot2 + theme(legend.position = "top",
                         panel.grid.major = element_blank(),
                         panel.grid.minor = element_blank(),
                         legend.title = element_text (size =9), 
                         legend.text = element_text (size=9), 
                         legend.key.width=unit(0.8,"cm")))
plot22

plot_cw1 <- grid.arrange(plot11, plot22, nrow = 1)
plot_cw1
ggsave(file="plot_cw1.png",dpi=600, width=10, height=5,plot_cw1)

describe(eu28$settlement_c)



plot3 <- interact_plot(model = cwmodel3, pred  = fcv_c,
                       modx  = settlement_c,
                       x.label = "FCV-19S",
                       modx.values = c(-0.29, 0.29),
                       modx.labels = c("Smaller","Bigger"),
                       y.label = "Climate change concerns",
                       legend.main = "Settlement size",
                       interval = TRUE)
plot33 <- (plot3 + theme(legend.position = "top",
                         panel.grid.major = element_blank(),
                         panel.grid.minor = element_blank(),
                         legend.title = element_text (size =9), 
                         legend.text = element_text (size=9), 
                         legend.key.width=unit(0.8,"cm")))
plot33

ggsave(file="plot_cw2.png",dpi=600, width=5, height=5,plot33)


#plotting random effects for countries#--------------------------------------------
  

#extracting per country slopes for COVID-19 concerns#--------------------------------------------

  
cntry_labels <- c("1" = "Austria", "2" = "Belgium", "3" = "Bulgaria", "4" = "Croatia", "5"="Cyprus",
                    "6" = "Czech Rep.", "7" = "Denmark", "8"= "Estonia", "9"="Finland", "10" = "France",
                    "11" = "Germany","12" = "Greece", "13"="Hungary", "14"="Ireland", "15"="Italy", "16"="Latvia", "17"="Lithuania", 
                    "18"="Luxembourg", "19"="Malta", "20"="Netherlands", "21"="Poland", "22"="Portugal", "23"="Romania", 
                    "24"="Slovakia", "25"="Slovenia", "26"="Spain", "27"="Sweden", "28"="United Kingdom")
  
  
##Per-country plots climate change concerns#------------------  

caterp_label <- c("(Intercept)"="Intercept","fcv_c"="Slope") ##labels
str(cw_ran <- ranef(cwmodel3)) ##extracting intercepts and slopes
cw_ran
str(cw_ran2 <- as.data.frame(cw_ran))
cw_ran2 


cw_ranef <- ggplot(cw_ran2, aes(y=grp,x=condval)) + xlab("FCV-19S") +
  ylab("Climate change concerns") + geom_point() + facet_wrap(~term,scales="free_x", labeller=labeller(term=caterp_label)) +
  scale_y_discrete(labels = cntry_labels) +
  geom_vline(aes(xintercept=0), colour = gray(1/2), lty = 2, size = 0.2) +
  theme_classic() +
  theme(strip.background = element_blank()) +
  geom_errorbarh(aes(xmin=condval -2*condsd,xmax=condval +2*condsd), height=0)
cw_ranef

ggsave(file="cw_ranef.png",width = 6, height = 5, dpi=600, cw_ranef)


##Per country coefs#------------------

cw_coef <- as.data.frame(coef(cwmodel3)$cntry)
cw_se <- as.data.frame(se.coef(cwmodel3)$cntry)
cw_merge <- bind_cols(cw_coef, cw_se) 
cw_merge$cntry <- c(1:28)  ##adding country number
cw_merge$cntry <- as.factor(cw_merge$cntry)
cw_merge <- within(cw_merge, cntry <- reorder(cntry, fcv_c...2)) 


cw_merge_mean <- cw_merge %>% ##extracting means
  summarize(mean_val = mean(fcv_c...2))
cw_merge_mean

cw_coef_plot <- cw_merge %>%
  ggplot(aes(y=cntry,x=fcv_c...2))+geom_point(aes(fcv_c...2,cntry))+xlab("FCV-19S coefficients") +
  ylab("Climate change concerns")+ scale_y_discrete(labels = cntry_labels)+
  geom_vline(data=cw_merge_mean, aes(xintercept=mean_val), colour = gray(1/2), lty = 2, size = 0.2) +
  theme_classic() +
  theme(strip.background = element_blank())+
  geom_errorbarh(aes(xmin=fcv_c...2 -2*fcv_c...21,xmax=fcv_c...2 +2*fcv_c...21), height=0)
cw_coef_plot


ggsave(file="cw_coef_plot.png",width = 6, height = 5, dpi=600, cw_coef_plot)


# percountry slopes in one plot

eu28$cw3 <- predict(cwmodel3)

cw_lines <- eu28 %>% 
  ggplot(aes(fcv_c, cw3)) + 
  geom_smooth(se = F, method = lm, size = 2) +
  stat_smooth(aes(color = cntry, group = cntry),
              geom = "line", alpha = 0.4, size = 1) +
  theme_bw() +
  guides(color = F) +
  labs(x = "FCV-19S", 
       y = "Climate change concerns", 
       color = "Country")+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),panel.border=element_blank())
cw_lines

ggsave(file="cw_lines.png",width = 6, height = 5, dpi=600, cw_lines)


##percountry slopes#------------------

eu28$cw3 <- predict(cwmodel3)

cw_slopes <- ggplot(eu28, aes(x = fcv_c, y=cw3)) +
  stat_smooth(method = "lm", colour="grey", level = .95) + 
  labs(y = "Climate change concerns", x = "FCV-19S")+ 
  facet_wrap(~cntry, labeller = labeller(cntry = cntry_labels),ncol = 4) + 
   theme(text=element_text(size = 7), axis.text.y = element_text(size=4), 
         axis.text.x = element_text(size=4),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank())
cw_slopes

ggsave(file="cw_slopes.png",width = 6, height = 5, dpi=600, cw_slopes)
 

